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ABSTRACT 

We propose an analytical 3-D model of the open field-line region of a neutron 
star (NS) magnetosphere. We construct an explicit analytic solution for arbi- 
trary obliquity (angle between the rotation and magnetic axes) incorporating 
the effects of magnetospheric rotation, relativistic flow of charges (e.g. primary 
electron beam) along the open field lines, and E x B-drift of these charges. Our 
solution employs the space-charge-limited longitudinal current calculated in the 
electro dynamic model of Muslimov & Tsygan (1992) and is valid up to very high 
altitudes nearly approaching the light cylinder. We assume that in the inner- 
most magnetosphere, the NS magnetic field can be well represented by a static 
magnetic dipole configuration. At high altitudes the open magnetic field lines 
significantly deviate from those of a static dipole and tend to focus into a cylin- 
drical bundle, swept back in the direction opposite to the rotation, and also bent 
towards the rotational equator. We briefly discuss some implications of our study 
to spin-powered pulsars. 

Subject headings: theory — pulsars: general 



1. INTRODUCTION 

The growing observational data (from radio- to 7— rays) on spectra and pulse pro- 
files of spin-powered pulsars prompt continued improvement of theoretical models (see e.g. 
Harding 2005, Kaspi et al. 2004). For example, study of particle acceleration and radiation 
produced at high altitudes in a pulsar magnetosphere (Muslimov & Harding 2004a [MH04], 
Hirotani et al. 2003) depends heavily on our knowledge of the structure of open magnetic 
flux lines (passing through the light cylinder) at very high altitudes, where the standard 
static magnetic dipole approximation is no longer accurate. Thus, the calculation reported 
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here was undertaken initially with the specific purpose of a quantitative description of the 
distorted magnetic field in a realistic magnetosphere at high altitudes . In this paper we 
develop the corresponding analytic solution which can be used up to very high altitudes 
nearly approaching the light cylinder. 

Since the pioneering theoretical studies (Goldreich & Julian 1969; Ostriker & Gunn 
1969; Sturrock 1971; Mestel 1971; Ruderman & Sutherland 1975; and Arons & Scharlemann 
1979) of pulsar phenomena there continues to be interest in the magnetospheric structure 
of a rotating neutron star (NS). The equation governing the structure of an axisymmetric 
pulsar magnetosphere was derived (under quite strong assumptions and idealizations) more 
than three decades ago (see e.g. Mestel 1973; Scharlemann & Wagoner 1973; Michel 1973; 
Okamoto 1974; Mestel et al. 1979 and references therein) and can be reduced to the well- 
known (special-relativistic) force-free Grad-Shafranov equation (see Grad 1967; Shafranov 
1966 for generic version of the equation). Since then the basic ideas of these and similar 
studies have shaped the school of thought that seeks to construct a mathematically closed 
(albeit highly idealized and axisymmetric) model of a NS magnetosphere and wind zone. The 
contemporary development of this school is mostly represented by Mestel and collaborators 
(see Mestel 1999 for general overview; and also Goodwin et al. 2004 for the most recent 
version) and by the Lebedev Institute group (Beskin et al. 1983). Recently, Bogovalov (2001) 
studied the MHD plasma flow in the magnetosphere of an oblique rotator with an initially 
split-monopole magnetic field. However, his solution cannot directly apply to radio pulsars 
since it is valid when R a < _R lc (where R a and i? lc are the Alfven and light cylinder radius, 
respectively). Thus, despite significant progress in the numerical solution (see Contopoulos et 
al. 1999; Mestel 1999) of the Grad-Shafranov pulsar equation, numerical simulation of plasma 
in a rotating NS magnetosphere (Biltzinger & Thielheim 2004; and Spitkovsky 2004), and 
also in the computation of MHD winds (Bogovalov 2001 and references therein), it is difficult 
to find in the literature any useful estimate of the high-altitude distortion of open-field lines 
of e.g. initially dipolar magnetic structure. How does the magnetospheric distortion at high 
altitude depend on the pulsar obliquity? Is the space-charge-limited current sufficient to 
distort open-field lines at high altitudes, and if so, how will it change the form of the open- 
field-line bundle? What is the high-altitude radial dependence of B? None of the existing 
NS magnetosphere models can readily provide clear and simple answers to these and similar 
questions. The main reason is that these models do not have simple analytic versions. The 
only available analytic model is the classical "vacuum" model of Deutsch (1955) which is 
not applicable (see also Section 3.1) to the physical situation in a real pulsar magnetosphere 
filled with charges and currents. 

In this paper we approach the problem in a slightly different way by first identifying 
and understanding the main physical effects distorting the geometry of the open-field line 
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configuration within the light cylinder. Then, taking advantage of the fact that these effects 
enter Maxwell's equations either as the first or second order terms in the radial distance, 
scaled by the light cylinder radius, we can significantly facilitate the problem by separating 
the multiple terms in the coupled system of equations. In doing so, we are able to solve 
analytically the simplified Maxwell's equations to determine the corrections to the static 
magnetic field caused by each of these effects. The main element of our model is our use 
of the electric current along the open field lines, which is determined by the space-charge- 
limited flow solution in the electrodynamic model of Muslimov & Tsygan (1992 [MT92]). We 
try to keep our formalism as simple as possible, so that our solution can be easily reproduced. 
The solution presented in this paper, being mathematically more transparent, allows us to 
understand the general picture of how the open magnetic field lines of a NS get distorted at 
high altitudes, depending on the pulsar obliquity, on the distribution of the electron current 
over the polar cap (PC) and therefore on the magnetic latitude and azimuth, and on the 
E x B-drift of charges. More importantly, the presented solution illustrates the way each 
of the effects of rotation, charge flow and E x B -drift contributes to the resulting pattern 
of distorted open field lines. Thus, the treatment is an attempt to explore a more realistic 
situation that is intermediate to the extreme cases of vacuum and MHD studied in the past. 
Although in the present study we focus on the region of pulsar magnetosphere confined 
by the light cylinder, we understand that the particles streaming along the open field lines 
will form a relativistic wind zone (Mestel et al. 1979). The regime of relativistic wind and 
corresponding global configuration of the magnetic field will be discussed separately. 

The paper is organized as follows. In § 2 we present a set of fundamental electrodynamic 
equations that will be employed in our study. In § 3 we formulate our approach and discuss 
how to incorporate the effect of rotation and relativistic charge flow (§ 3.1), and the effect 
of E x B-drift of electrons (positrons) on the structure of open field line region within the 
light cylinder of a NS magnetosphere (§ 3.2). In § 4 we provide 3-D views of open field lines 
for different obliquities. Finally, in § 5 we discuss the results of our study and summarize 
its most exciting implications for spin-powered pulsars. 



2. Basic Equations 

Let us consider the magnetosphere of a rotating NS and assume that in the frame of 
reference rigidly corotating with the NS the magnetic field is stationary. The very gen- 
eral equations describing the electromagnetic field produced by the rotating NS in the Lab 
(inertial) frame are the first couple, 

V-B = 0, (1) 
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VxE = — — , 

c at 

jquations, 



(2) 



V • E = 4irp, 



(3) 



V x B 



1 <9E 4tt . 

c dt c 



(4) 



where p and j are the electric charge and current density, respectively. 

In this paper we adopt the standard picture that the NS's magnetosphere has two 
distinctive regions: a 'dead zone' with field lines that close within the light cylinder, and 
without any current flow along the field; and the open-field line region extending beyond 
the maximum radius of corotation (in cylindrical coordinates with the z-axis along the NS's 
rotation axis), the light-cylinder radius R\ c c/fi, and with flow of charges along the field 
lines. As in our previous studies, we will be working in a spherical polar coordinate system, 
[r] = r/R, 6, 0), in which the polar axis is parallel to the magnetic moment. In this system, 
r is the radial coordinate and R is the stellar radius, 9 is the polar angle measured from 
the magnetic dipole axis, and is the azimuthal angle measured counter-clockwise from the 
meridian passing through the rotation axis. We refer to this coordinate system throughout 
the paper as magnetic coordinates. Finally, we define x to be the pulsar obliquity (angle 
between the NS rotation axis and magnetic dipole moment). 

In our model calculation we will assume that the static (unperturbed by rotation and 
currents) magnetic configuration of the NS has a pure dipole geometry, 



where e r , e e are the corresponding basis vectors of the magnetic coordinate system defined 
above and is the magnetic field strength at the magnetic pole. Here, for the sake of 
simplicity, we ignore static general relativistic corrections to the magnetic field. This is 
well justified, because we are interested in the corrections to the magnetic field at high 
altitudes approaching the light cylinder and caused by the magnetosphere rotation and 
flow of charges along the open field lines. The only exception will be the expression for 
the longitudinal component of the current density (see e.g. formula [27] below) which is 
essentially determined by the condition at the stellar surface where general relativistic effects 
are not merely important but make a qualitative difference. 




(5) 
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3. Formulation of the problem and analytic solution 

3.1. The effect of rotation and charge flow 

We will be searching for the steady-state solution to equations (l)-(4). In this case, 
the time derivatives in equations (2) and (4) are determined by the rotation of the NS 
magnetosphere relative to the Lab frame. It is important that well within the light cylinder 
(i]~Vic = R\c/R) this rotation is most likely a solid-body rotation, and we can use the 
following transformations of partial time derivatives between the Lab frame (subscript "Lab" ) 
and the frame of reference corotating with NS magnetosphere (subscript "corot"): 

dB\ (dB^ 



J corot V J Lab 



(dE) (dE) „ . ^ x „ „ 



(7) 

at l -\dt\ +u "*- v " ( 8 ) 

Ul ) corot l UL J Lab 

where u rot is the rotational velocity of the magnetosphere. 

We assume that in a steady state, the time derivatives in the LHS of equations (6)- (8) 
vanish, so that Maxwell's equations (2) and (4) can be rewritten in the following form 

VxE = -Vx(/3 rot xB), (9) 

An 

V x B = V x (/3 rot x E) - /3 rot V • E + — j. (10) 

where (3 Iot = u rot /c. 

Also, the charge continuity equation, 

! + V.j = 0, (H) 
with the help of relationship (8) takes the form 

urot- Vp- V-j = 0. (12) 

By combining equations (3) and (10), we get 

V x(B-^xE) = -(j- pu rot ). (13) 

V c / c 
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Thus, the steady-state solution for a rotating magnetosphere is determined by equations (1), 
(9), (12) and (13). Note that, since V • u rot = and V • (pu rot ) = u rot • Vp, equation (12) 
translates into 

V-(j-pu rot ) = 0. (14) 

Now we should discuss the physical origin of the current density j. Within the light cylinder 
the current is mostly determined by the longitudinal (owing to the relativistic electrons 
streaming along the magnetic field lines) and rotational (owing to the bulk rotational motion 
of charges) components, 

j=j||+jrot- (15) 

In this Section we ignore the effect of E x B-drift. This effect will be discussed separately 
in Section 3.2. The main reason is that the net current produced by the E x B-drift may or 
may not vanish depending on the specific scenario of particle acceleration within the open 
field line region. For example, if the relativistic beam is quasineutral, then charges of both 
sign will be drifting in the same direction and with the same velocity thus producing zero net 
current. On the contrary, if the beam is charged (e.g. primary electrons and quasi-neutral 
electron-positron plasma are flowing in the region with open field lines), then the E x B-drift 
of electrons (positrons) can significantly contribute to the net current and therefore affect 
the structure of the magnetic field at high altitudes (see Section 2.2 for details). 

The longitudinal and rotational components of the current density can be written as 

Jii=in§ (is) 

and 

jrot = -|p|u rot , (17) 

respectively. 

Note that in equation (17) the "-" sign signifies that negative charges (electrons) are 
involved in rotational motion. 

By using expressions (16) and (17) we can rewrite equation (13) as 

Vx(B-/3 rot xE) = ^j,|. (18) 

Although this equation is essentially the same as the corresponding equation derived by 
Beskin et al. (1983), we should point out that there are some principal differences: Beskin 
et al. neglected the component j rot in their expression for j; and, in their derivation of the 
equation (see their eq. [15]) analogous to our eq. (13), they neglected the component pu rot 
resulting from the transformation (7). 
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We assume that at higher altitudes where the effect of rotation becomes increasingly 
important, the electric field (in the Lab frame) is mostly determined by rotation (see also 
MH04), 

E » -/3 rot x B d . (19) 

Note that, besides j ^ and p ^ 0, relationship (19) assures the fundamental difference 
between ours and Deutsch's solutions. Also, we should point out that for high altitudes, 
the classical vacuum analytic solution of Deutsch transforms into a pure wave-like solution 
well beyond the light cylinder and is hardly applicable to the realistic situation. Finally, 
the Deutsch's solution is presented in spherical coordinates with the polar axis along the 
rotation axis (see also Cheng et al. 2000), whereas our solution will be presented in magnetic 
coordinates. 

Then, equation (18) can be rewritten as 

VxB = — j,|+J, (20) 

c 

where 

J = V x G, (21) 

and 

G = (3 TOt xEe fi ot B d - (3 TOt ((3 TOt • B d ). (22) 
Here we assumed that E is determined by equation (19). 

The main goal of the present study is to construct the appropriate analytic solution for 
the open magnetic field lines valid within the light cylinder. By inspecting the Maxwell's 
equations derived above one can see that the corresponding solution for vector B can be 
generally presented as 

B = B d + B (1) + B (2) + B (2 *) , (23) 

where B d is a pure dipole component (see formula [5]), B^ is the first-order correction to 
the static dipole component which is ~ (v/Vic) B d , and B^ 2 ) and B^ 2 *) are the second-order 
corrections, ~ (v/Vic) 2 B d , respectively. The correction B^ is produced by the charge flow 
along the open field lines. B^ 2 ) is the distortion caused by rotation of B d , and B^ 2 *) is 
the distortion generated by E x B drift of the outflowing charges. The terms B d and B^ 2 ^ 
represent a rotating vacuum solution subject to the condition of equation (19). 

In what follows, for the sake of convenience, we will use the magnetic spherical coordi- 
nates with the polar axis along the magnetic dipole moment. 

Obviously, the first-order correction to the dipole magnetic field is generated by the 
longitudinal current flowing along the poloidal (and mostly determined by the dipolar com- 
ponent) magnetic field. In this case the contribution from the displacement current J is of 
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the second order and can be neglected. Thus, the equation for determining reduces to 

4-7T 

V x B (1) = — jii, (24) 

c 

To complete the formulation, we should add the following couple of equations (see equations 
[1] and [14]) 

V • B (1) = 0, (25) 

and 

V-j|,=0. (26) 

To solve the system of equations (24) - (26), we need the explicit expression for jy. For this 
purpose we will employ the electrodynamic model of MT92 and write 

,B d n 



JII = ~c\p\ 



B d ~ 2tt 



3 

(1 - re) cos x + ^ d o f sin X cos 



B d , (27) 



where 9$ ~ (VtR/c) 1 / 2 is the canonical PC half-angle, £ is the dimensionless magnetic colati- 
tude of open field lines (£ = 1 corresponds to the last open field lines, and £ = corresponds 
to the magnetic axis), and re is the parameter measuring the general-relativistic effect of 
frame dragging at the stellar surface in units of stellar angular velocity Q. According to 
our estimate, for most more or less realistic NS equations of state, re « 0.15 where 
I45 = J/10 45 g-cm 2 , R 6 = R/10 6 cm, I is the moment of inertia of NS of radius R. 

One can easily verify that expression (27) for jy satisfies equation (26). Note also that in 
formula (27) the explicit ^-dependence implies that for pure dipole field lines the azimuthal 
coordinate coincides with the azimuthal coordinate of a stream line. In other words, for 
any given value of one can calculate the current density jy which is fixed at the stellar 
surface at the azimuth 0p C (= 0). 

Let us now introduce the dimensionless vector, b/ 1 ), such that 

B(1) =(v) 5o d b«(^,0), (28) 
where Bq is the dipole field strength at the magnetic pole. Then, equation (24) reduces to 
^{bf ^8)-^bP = -l(a + Pcos<j>)smecose, (29) 



09 \ * J d(f) " rf 

-J^^ 6 ?) =-^(a + /3cos0)sin£, (30) 

0. (31) 



1 i)h}' 
sin 6 d(f) 

d ( ,nA db^ 



drjV " J 09 
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where a = (1 — k) cosx, and (3 = (3/2)# £ sin \- 

By examining the system of equations (29) - (31) and (25) one can see that the solution 
for that vanishes at infinity, should have the following simple form 

^ = ^sin0, (32) 

,(i) g{9) . , /qq\ 

% = — 5- sm0, (33) 



b { l ] = ^i»i(8) + h 2 (9) cos 0]. (34) 

This solution should be regular at the magnetic pole (at 9 = 0) and satisfy the periodic 
condition, b^(0) = b (1 )(0 + 27r). By substituting expressions (32)-(34) into equations (25), 
(29)-(31), we find that equation (29) is a consequence of equations (25), (30) and (31). Also, 
the solution for hi{9) can be written immediately, 

h^O) = -a sin9, (35) 

whereas solution for /, g, and h 2 can be found from the system 

/" sin 2 9 + f sin 9 cos 9 - f = (3 sin 2 9, (36) 

9 = -/', (37) 

h 2 = -J—-Psm9, (38) 
sin 9 

where ' = d/d6. The solution of equation (36) that is finite at 9 — > reads 

/ = (3(1 -9 cot 9), (39) 

so that g and h 2 can be easily determined after inserting this solution into equations (37) 
and (38), respectively. 

Thus, the final analytic solution for can be written as 

B d 

| # f (1 -#cot#)sinxsin0, (40) 
77^ £?o f 9 — sin 9 cos N 



1 



M 1} 









n e — 


3 
~2 




\ B$_ 


(1- 




J ^3 



I rt S \ -2/1 

i]\ c J \ sin 

„ 3„ /sin 3 + sin0-0 cos0 



sin x sin 0, (41) 



2 U M sin 2 



sin x cos 



(42) 
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where 




and 



A = (1 - k) cos X + (3/2) 9 £ sin X cos0°. (45) 



Note that A (cf. expression [27]) depends on 0p C not 0, simply because formula (44) implies 
that the stream lines of current determining B^ 2 ) now have an azimuthal component, in 
which case the coordinate <fi ^ 0p C and can only be set by specifying its value (or the 
value of A ) at magnetic azimuth 0° c at the PC surface. Here, for the sake of simplicity, we 
assume that in formula (45) the parameters are such that A > 0, i.e. the charges of the 
same sign (electrons) can be ejected from the stellar surface. The case of Ao < is briefly 
discussed in the last paragraph of Section 4. 

Obviously, we can search for a solution for the component B^ 2 ) as having the following 
dependence of the dimensionless vector h^ 2 \ 



where the spherical components of vector g are given in Appendix A. 

By using expressions (44), (47) and employing the solution for B^, equation (43) can 
be reduced to 




(46) 



Note also, that the vector J can be presented as 




(47) 



V x b (2) = V x (g + h), 



(48) 



where the spherical components of vector h are also given in Appendix A. 
From equation (48) we can get 

bf ] =g r + h r . 

The 9— and 0— components of b( 2 ) can be presented as 



(49) 




(50) 
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l(2) 



g<p + he/, + 



l OX C 2 (0) 



+ 



(51) 



77 sin 6* (90 77 sin # 

Thus, b^ 2 ) vanishes at large radial distances, and X(9, 0), Ci(6>), and C 2 (<fi) are some functions 
to be determined from equation 

V • b( 2 ) = 0, 



(52) 



which translates into 



r}(g r + h r ) + 



d 



sin 9 1 89 
1 



sm9 [ r]g e + Ci + — J 



+ 



sin6> 



d , . i d 2 x dc 2 



0. 



(53) 



By separating the variables in equation (53) and assuming the regularity of solution at the 
magnetic pole and periodicity over azimuthal coordinate (see also the comment following 
equation [34]), we can find that 

Ci = -jjj [sin 2 x (1 + cos 2 6») + 2 cos 2 x sin 2 6»] + A (1 - k) cosxj sintf, (54) 



and 
where 



C 2 = 0, 

X — (U sin x cos x + V A $o £ sin x) cos <fi + W sin 2 x cos 20, 



U 



— !— (cos 9-1)-- sin 9 (cos 2 0-^1, 
sin# v ; 2 V 2 ' 



f cos0 + ^4- -4 
2 1 sin# 



and 



W = - 
8 



5 /l-cos0 1 



sin 2 9 cos 9 



.2 V sin 2 2, 
By using the above expressions, we can now write the components of B^ 2 \ 

B[ 2) = (jj-^J ^| cos9 |cosx[cosx sin 2 9 + 2 A (1 - «)] + sin 2 x ^1 - ^ sin 2 ^ 



(55) 
(56) 

(57) 
(58) 

(59) 



3 (9 

cos x sin 6> cos + - A 6> f I -— 

2 V sm cos 9 



sinx cos0- 



- sin 2 x sin 2 cos 20 \ , 



(60) 
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B 



(2) 



vY B* . 
— — smU 

Vic J V 



— cos x sin 2 & + Aq (1 — «) 



cosx- 



1 sin 2 X (l " ^ sin 2 ^ + 

3 f 3(l-fl cotfl) 

2 V sin p 



1 - COS 9 1 \ l.o. 

sin^V^T^ + 5 Sm ' 



1 n 1-COS^^ 

- cot 6^ cosztH o I cosx- 

4 sin 6* 



sinx cos0— 

sin 2 x cos 20 )■ , 



(61) 



5 



(2) 



rj_ 

.Vic 



2 R d 

rf 



1 l-cos0\ 3 A „ , 3 sin# cos9-9 (3-2 sin 2 #V 

7+ -2, cosx + -A ^e : Yn ' 

4 sm 9 1 2 sm 9 



x 



sm x sm — 



5 / 1 — cos 9 



1 



8 V sin J # 2 sin# 



sin x sm 20 f , 



(62) 



3.2. Effect of E x B drift 

The current density associated with the electron drift within the region of open field 
lines can be written as 

j ExB = -|Po|^ExB, (63) 



JexB 

where 



l£2" 



IPol = ^- Ao 4> (64) 

2llC V 

from equation (27). By using the approximation (see [19]) E — /3 rot x B, we can rewrite 
expression (63) as, 

JexB = -|Po| C /^rot + |Po| C /3 ro t,|h ( 65 ) 

where /3 rotj|| = (B d /5 2 ) (/3 rot • B d ), and 5 » 5 d 

The second term in (65) adds to the component of current density jy oc B d (see eq. [27]). 
Here we have a situation where non-relativistic (or even nearly relativistic) drift motion in the 
longitudinal direction (along the magnetic field lines) is superposed on essentially relativistic 
flow. Thus, without any loss of generality, we can justifiably assume that this component 
of drift motion can be ignored in the longitudinal relativistic flow of charges. However, the 
first term in j ExB should be explicitly added to = —\po\ c B^ / B (see equation [43]), so 
that the corresponding correction to the magnetic field, B^ 2 *\ produced by the E x B-drift 
current will be determined by the equation 

VxB( 2> > = ^f», (66) 
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where 

j( 2 *) = -| Po |c/3 rot (67) 
is the current density associated with the E x B drift of electrons. 

Now we can introduce the dimensionless vector b>( 2 *) via the formula 

B (2* ) = S d Aob (2*) (68) 

Then equation (66) translates into 

R V x b (2 * } = i (2 *\ (69) 
where vector i^ 2 *) has the components 

if* ] = 0, (70) 



2 



H = -o sin X shi0, (71) 

and 2 

^(2*) _ — (cos% sin 9 — sin x cos 9 cos0). (72) 
i] 2 



The analytic solution for equation (69) satisfying V • t>( 2 *) = can be easily written as 

2 

^(2*) _ — (cosx cos 9 + sin x sin# cos</>), (73) 
bf*" 1 = - (cos x sin 9 — sin x cos 9 cos 4>) (74) 

and 

6*?*^ = - sin x sin 0. (75) 

Thus, the components of vector B^ 2 *) can be explicitly presented as 

(n \ 2 B d 
— -4 A (cosx cos # + sin v sin 9 cos0), (76) 
VicJ V 3 



2 Sn d 



i?! 2 *- 1 = ( — ) —j- A (cosx sin ^ — sin x cos# cos0) (77) 
VW V 



and 



2 5n d 



B?* ) = (- ^A o sin X sin0. (78) 
VW T 

Now, by inserting the corresponding components determined by formulae (5), (40) - (42), 
(60) - (62), and (76) - (78) into the general expression (23), we can calculate the structure of 
the open magnetic field lines all the way from the NS surface up to the high altitudes nearly 
approaching the light cylinder radius. 
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4. Visualization of Analytical Formulae 

By using the spherical components of B given by the general expression (23), we can 
easily produce a 3-D plot of any open field line for arbitrary obliquity In Figures 1-3 we 
depict some views of a set of open field lines emanating from the same magnetic colatitude 
of the NS PC and equally spaced (by 30°) in magnetic azimuth. We use the Cartesian 
coordinates X, Y and Z with the center at the magnetic pole. Also, the Z-axis is along the 
magnetic moment, and the positive X-axis is pointing towards the rotation axis. Finally, all 
the axes are scaled by the light cylinder radius, and a pulsar spin period of 0.3 s is used. In 
Fig. 1 we show the side views of open field lines having dimensionless magnetic colatitude 
£ = 0.9 for different obliquities x — 0> 30, 60, and 90°. [Here, for the sake of simplicity, we 
normalize the magnetic colatitude by the canonical PC half-angle, 9q. In general, the foot 
points of last open-field lines should be determined by tracing back from the light cylinder 
to the stellar surface along the last open-field lines. For the purpose of the present study this 
exercise is not essential, though.] In Fig. 2 we present the same set of magnetic field lines 
as in Fig. 1 but viewed from the magnetic pole. For the aligned case (x — 0), as one can see 
from Figures 1 and 2, the field lines are axisymmetric and get swept in the direction opposite 
to rotation, in contrast to the vacuum case. The non-vanishing component (see equation 
[81] below) gives a sweep back due solely to the real current, even in the aligned case. For the 
aligned rotator the expressions (40) - (42), (60) - (62), and (76) - (78) significantly simplify, 
so that for the components of B defined by formula (23) we get 

B d 

B v = ^- cosfl{l + (77/77 lc ) 2 [sin 2 fl - 2 k (1 - k)]}, (79) 
B e = m sin#{l - (1/2) (r^ic) 2 [sin 2 9 An (1 - «)]}, (80) 

2 T]^ 

B d 

B^ = £ (l-«) ( V /rn c ) sin0. (81) 

From formulae (79), (80) one can see that at 9 < 9 h (where 6 h « [2«(1 - k)] 1/2 « 30°) the 
open field lines have slightly more flaring than those of a pure dipole, whereas at 9 > 9 h 
they are more focused towards the magnetic axis. Also, in a small-angle approximation 
(9 <C 1), from equations (79), (81) we can write the following approximate field-line formula, 
~ 0p C — (1 — k) (v/rfrc)- 111 Figure 3 we illustrate the effect of sweep-back for the aligned 
rotator and for the field lines emanating from different magnetic colatitudes (£ = 0.2, 0.4, 0.5, 
and 0.9). This Figure shows that the effect of sweep-back vanishes towards the magnetic 
axis. Note that if terms depending on the current (last terms in equations [79] and [80] 
and the whole RHS of equation [81]) are turned off, there still remains a contribution from 
displacement current. This contrasts to the Deutsch solution, where the displacement current 
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oc dE/dt oc <9E/<9(/> is zero in the aligned case. This is because the wavelike solution imposed 
by Deutsch at large distances requires non-axisymmetry in order to produce a not vanishing 
displacement current. 

As the obliquity increases (see cases x — 30 and 60° in Figures 1 and 2) the leading 
(negative Y) and trailing (positive Y) field lines become asymmetric (with respect to rotation 
by 180 ° around magnetic axis). Also, in this case (Fig. 1, x = 30 and 60°) and at high 
altitudes, the field lines get more focused and the entire bundle bends away from rotation 
axis. For the orthogonal rotator (case x = 90° in Figures 1 and 2) and for the value of 
spin period (0.3 s) we used in our numerical calculation, both the effect of relativistic flow 
and NS rotation significantly diminish, so that the configuration of open field lines (at least 
up to ~ 0.3 ?7i c ) is practically the same as in the case of a static magnetic dipole. This is 
difficult to illustrate for arbitrary values of 9 for open field lines, because for x — 90° the 
analytic formulae still look rather cumbersome. However, in a small-angle approximation, for 
X = 90°, it is easy to see that B T ~ Bf, Be ~ By, and ~ 9 (r)/r)\ c ) 2 Bf, which explains 
the patterns depicted in Figures 1 and 2 for the case of x — 90°. This behavior sharply 
contrasts to the case of the vacuum orthogonal rotator, where the field-line sweep-back is a 
maximum (Arendt & Eilek 1998; Cheng et al. 2000; and Dyks & Harding 2004). 

Note that the components of the magnetic field and B^ 2 *) depend on the "source" 
function Ao, which is determined by the obliquity, spin period (through 9o), magnetic co- 
latitude and azimuth of a field-line foot point at the PC surface. For electrons accelerating 
from the PC along the favorably curved (cos</>p C > 0) field lines, A is always positive (we 
assume that < x < 90°, i.e. "normal polarity" pulsar; see also MH04). However, for the 
unfavorably curved (cos0p C < 0) field lines there may be a situation where the second term 
in A dominates. In this case (e.g. for the millisecond pulsars with high obliquities) the un- 
favorably curved field lines may become inefficient in providing continuous steady-state flow 
of electrons (see also Muslimov & Harding 2003, MH04). This effect should be taken into 
account in calculating the 3-D magnetic structure for the short-period pulsars. We should 
also point out that one of the advantages of having our analytic expressions for and 
B( 2 *) depend on the source function A is that the latter can be replaced by e A , where 
< e < 1 is an arbitrary factor which takes into account the possibility that at the stellar 
surface the electron current may be a factor of ~ e less than the Goldreich- Julian value. 
This is important for modeling the effect of relativistic electron flow of different magnitude 
on the magnetospheric structure. In other words, by using our model we can explore the 
space-charge-limited flow approximation (e.g. by measuring the magnitude of the effect of 
field line focusing and sweep-back in pulsars), examine the role of E x B drift in determining 
the magnetospheric structure at very high altitudes, and probe the occurrence of different 
acceleration conditions on favorably (cos0p C > 0) and unfavorably (cos0p C < 0) curved field 
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lines. We plan to address these and other consequences of our model, including possible 
observational tests, in subsequent studies. 

5. Discussion and Conclusions 

In this study we began a quantitative analysis of the distortion of open magnetic field 
lines of a rotating NS with an arbitrary obliquity angle in the presence of relativistic charge 
flow. The static field configuration is assumed to be dipolar. Our analysis is aimed at the 
derivation of simple analytic formulae which can be used up to very high altitudes, say 
~ 0.5 — 0.7 of the light cylinder radius. We presented the explicit analytic expressions for 
the first and second order corrections to the static dipole magnetic field that are produced 
by the effect of relativistic flow of charged particles (e.g. with the negative net charge) 
along the open field lines, bulk magnetosphere rotation, and non-vanishing E x B-drift of 
a net charge across the open field lines. For the longitudinal current, we employed the 
corresponding expression derived earlier in the space-charge-limited-flow approximation by 
MT92. The longitudinal current produces a substantial distortion of the open field lines 
at high altitudes by twisting them in a direction opposite to that of rotation (as viewed 
down the magnetic pole) and making the entire bundle of open field lines less flaring (more 
focusing) along the magnetic axis (see equations [40]- [42]). The effect of field-line twisting 
is clearly seen for small obliquities (see e.g. Fig. 2, cases x = 0° and 30°), in which case it 
can also be recognized as the sweep-back effect. The effect of focusing of open field lines can 
be seen in Figure 1 (cases x — 30° and 60°). The magnitude of these effects is determined 
by the obliquity and spin-period, so that these effects are much more pronounced for small 
obliquities, in contrast to the vacuum case. The global current associated with the bulk 
magnetosphere rotation distorts the open field lines (see terms that are not proportional to 
A in equations [60]- [62]) in such a way as to bend (at very high altitudes) the bundle of 
open field lines down toward the rotational equator. Finally, the effect of E x B -drift of 
charges (of the same sign as the net charge of relativistic longitudinal flow) across the open 
field lines results in some asymmetry between the leading and trailing edges of the bundle 
of open field lines (see equations [76] -[78]). The distortion of open field lines caused by this 
effect can be seen in Figure 1 (cases x — 30° and 60°). 

A number of models of emission from rotation-powered pulsars have relied on the mag- 
netic field structure of the vacuum retarded dipole (Deutsch solution). The predicted high- 
energy pulse profiles in outer gap (e.g. Cheng et al. 2000), two-pole caustic (Dyks & Rudak 
2003, Dyks, Harding & Rudak 2004) and slot gap (MH04) models, where most of the emis- 
sion occurs in the outer magnetosphere, depend sensitively on the structure of the field at 



-17- 



high altitude. The large differences in the rotational distortions of the magnetic field in vac- 
uum and non- vacuum cases, that is demonstrated by our solutions, could produce important 
changes in the predictions of such models. 

Although the main area of application of our study was meant to be the modeling of 
particle acceleration and emission of energetic photons at high altitudes in the open-field line 
regions of pair-starved pulsars (especially in millisecond pulsars), the solution presented here 
should also be applicable to the majority of pulsars producing high pair multiplicity in their 
magnetospheres. This is because our solution assumes a primary current given by a space- 
charge limited flow model. Even though this model (and any acceleration model) necessarily 
departs from the force-free (ideal MHD) case, the departure is small. Furthermore, we have 
illustrated in previous studies (MH04, Muslimov & Harding 2004b) that even in the charge- 
starved limit, the actual space-charge along the open field lines approaches the Goldreich- 
Julian charge, and thus the force-free condition, at high-altitudes in the magnetosphere. 
Thus the difference between the high-altitude field structure of a pair-starved pulsar and a 
pulsar producing an abundance of pairs should be minimal. 

It is worth mentioning some other areas where our formulae can be used and subjected 
to observational tests. First of all, our analytic expressions can be implemented in the 
analysis of pulse polarization properties in radio pulsars to probe e.g. the geometry of the 
magnetic field lines in emission cites (e.g. Gangadhara & Gupta 2001; Blaskiewicz, Cordes & 
Wasserman 1991; Dyks, Rudak & Harding 2004; and Hibschman & Arons 2001). Second, the 
results of our study can be used for the interpretation of multi-frequency (e.g. in radio-, IR-, 
optical-, X-, and 7-rays) light curves of pulsars to explore the 3-D picture of pulsar emission 
(see Cheng et al. 2002; and Romani 2002 for a review). And third, our formulae might be 
useful in modeling of magnetosphere geometry, in general, and particle acceleration regions, 
in particular, that is needed in interpretation of observations. For example, the recently 
discovered double pulsar system PSR J0737-3039 A&B (Burgay et al. 2003; Lyne et al. 2004) 
may provide us with the opportunity to probe both the pulsar magnetospheric structure and 
wind density (see e.g. Arons et al. 2005; Rafikov & Goldreich 2005; and Zhang & Loeb 
2004). Also, our study may be applicable to the interpretation of pulsar breaking indices 
(see e.g. Manchester & Taylor 1977 for a general discussion). Finally, it is interesting to 
point out that in the course of pulsar evolution their spin periods and maybe field strengths 
(and/or inclination angles) as well as (as our present study implies) the geometries of their 
open-field line regions can change. This effect should be taken into account in population 
synthesis. 

In the future, we can extend our present study in, at least, three different ways: a) To 
derive the appropriate analytic solution valid beyond the light cylinder; b) To incorporate 
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the deviation of open-field line structure at high altitudes into the electrodynamic model 
to derive a more accurate high-altitude E\\; c) To proceed with detailed modeling of high- 
altitude pulsar emission and comparison with the observational data. Note that items b) and 
c) imply some additional development of our model, while item a) may turn into a separate 
study falling under the topic of a pulsar relativistic wind. In addition, after implementing 
item b) we will be able to perform well-founded validation of the existing scenarios of particle 
acceleration and emission in pulsars. Finally, as a result of effort c), we expect to come up 
with some constraints on any theoretical model of open field lines near the light cylinder. 

We acknowledge support from the NASA Astrophysics Theory Program through the 
Universities Space Research Association. 
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A. Components of vectors g and h 



In this Appendix we present explicit expressions for the spherical components of vectors 
g and h, which are defined in Section 2.1 (see formulae [47] and [48]). 

The r— , 9— and <fi— components of vector g: 

1 



and 



g T = - (cos x sin $ + sin X cos + sin y sin 9 sin — 
2 sinx cosx sin 9 cos 9 cos0) cos#, 

go = — (cosx sin 6 — sin \ cos 9 cos</>) 2 sin 9, 
zrj 

gs = — (cosx sin 6 1 — sin x cos 9 cos0) sinx sin 9 sin (p. 
2r) 



respectively. 

The r— , 9— and 0— components of vector h: 



h T = - A 
77 



'1 - re) cosx cos6» - - 9 £ 
h e = 0, 



sm 



9 cos 9 



sin 9 



■ sm x cos < 



and 



3 . „ _ 2 cos 9 + 9 sin - 2 

^ = - A 6*0 4 ^-7 sm x sm < 

i] sm9 



(Al) 
(A2) 

(A3) 



(A4) 
(A5) 

(A6) 
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Fig. 1. — The side views of a bundle of open field lines with £ = 0.9, for different values 
of pulsar obliquity \ ( = 0, 30, 60 and 90°). The Cartesian coordinate system (X, Y, Z) 
is centered at the magnetic pole, with the Z-axis along the magnetic moment, and positive 
X-axis pointing to the rotation axis. The coordinate values are in units of the light cylinder 
radius. The calculations are performed for the pulsar spin period of 0.3 s. Black lines denote 
favorably curved field lines (cos 0p C > 0) and gray lines denote unfavorably curved field lines 
(cos0° c <O). 




Fig. 2. — The views down the magnetic pole (from far above) of a bundle of open field lines 
with £ = 0.9 for different values of pulsar obliquity x (— 0, 30, 60 and 90°). Same as in 
Figure 1. Leading field lines are at negative Y and trailing field lines are at positive Y. 




Fig. 3. — The side views of bundles of open field lines with £ = 0.2, 0.4, 0.5, and 0.9 for the 
aligned case (x = 0). Same as in Figure 1. 



